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; H Floating Point Square Root routine 16-SEP-1984 01:40:3 AX/VMS Macro V04-00 Page 1 
. Hits 94:59:2 MTHRTL.SRCJMTHHSQRT.MAR; 1 . (1) 
TITLE MTHSHSQRT 4 floating Point Square Root routine 
: IDENT /1-001/ : File: MTHHSQRT.MAR EDIT: JAW1001 


MARAARAAAASAALALALALAALE LALLA LATTE ESSER EEE ERE EEE EERE PRESET TE TTT Tea 


COPYRIGHT (c) 1978, 1980, 1982, 1984 BY 
DIGITAL EQUIPMENT CORPORATION, MAYNARD, MASSACHUSETTS. 
ALL RIGHTS RESERVED. 


THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED 
NLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE 
St tte et THE ABOVE Sg lth yt f THIS SOFTWARE OR ANY OTHER 


SSSSssssssssss 


SOoooooooooooooo 


R 
TRANSFERRED. 


THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE 
PORPORAT ibn NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT 


DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS 
SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL. 
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; FACILITY: MATK LIBRARY 


:++ 
; ABSTRACT: 


; MTHSHSQRT is a function which returns the H-floating square root of 
; its H-floating argument. The call is standard call-by-reference. 


; MTHSHSQRT_RB is a special routine which is the same as MTHSHSQRT 
3 except that a faster non-standard JSB call is used with the argument 
in RO/R3 and no registers are saved. 


VERSION: 1 
> HISTORY: 
> AUTHOR: 
; John A. Wheeler, 25-Jul-79: Version 1 


MODIFIED BY: 
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TH ‘ fin Tbk iis Point Sayere Roos ene! im one 91: $9: i AX/VMS ee ro V04- On? mars Page 


Deta wrrent t dit H + MTHRTL. SRE JMTHHSQ 
§ -SBTTL HISTORY 3; Detailed Current Edit History 
3 Edit History for Version 1 of MTHS$HSQRT 
0 
1 


1-001 - Adapted from MTH$SQRT version 1-013 praesent version 1-013, 
and MTH$GSQRT version 1-001. JAW 25-Jul-79. 
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; H Floatin Point uare Root ro a 16-SEP- 1:40:3 AX/VMS Macro V04-00 P 
BECLARATIONS ; Dec erat ve Part ° Modu 6-SEP- 138 1:25: rk} YATART SRC JMTHHSQRT.MAR; 1 a 3) 
: ~SBTTL DECLARATIONS : Declarative Part of Module 
66 ; INCLUDE FILES: 
6 ; EXTERNAL SYMBOLS: 
29 7 DSA at. GBL ; Declare all externals 
7 “EXTRN MTHSSSIGNA 3 SIGNAL SEVERE error 
8 if -EXTRN MTHSK » SQUROONEG Error code 
74; 
° 75 ; EQUATED SYMBOLS: 
0 6 : 
000 7 
000041FC 0000 78 ACMASK = “M<IV, R2, R3, R4, RS, R6 R7, R8&> 
000 79 | hoatéeer save mask and IV enable 
Bae 80 
0000 81; 
000 Ht 3: MACROS: None 
B08 3 3 
O60 ae 3; PSECT DECLARATIONS: 
00000000 36 -PSECT _MTHSCODE PIC,SHR,LONG,EXE,NOWRT 
0000 8 ; Program section for math routines 
0000 88 
0000 89 ; 
0000 90 ; OWN STORAGE: 
0000 91; 
0000 3 3 
444 82 3; Constants A and B chosen for k odd 
13¢D4054 0000 95° LF_ODD_A = “X13CD4054 : Decimal: 0.8284271 
3C4A3F 98 O88 38 LF-ODD-B = “X3C4A3F98 : Decimal: 0.2973349 
444 4 : Constants A and B chosen for k even 
F61A4015 0000 100° LF_EVEN_A = “XF61A4015 : Decimal: 0.5857865 
4B8233FD7 0000 #£101 LF “EVEN, -B = “X4B233FD7 : Decimal: 0.4204951 
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1-001 : 6-SE gf o133 a: f ¢ 


Float 04- p 
MTHSHSQRT = Standard H-Floating SQRT  6=SEP-19 00 age 4 


V 
MTHRTL.SRC JMTHHSQRT.MAR; 1 (4) 
-SBTTL MTHSHSQRT = Standard H-Floating SQRT 


{FUNCTIONAL DESCRIPTION: 

3 HSQRT = H-floating point square root function 

; HSQRT(X) is computed using the following approximation technique: 
If X < 0, signal error. If X = 0, return HSQRT(X) = 0. 

Let X = 2e*K * F where F is the fractional part. 


If K is even, X = 2%*(2P) * F 
HSQRT(X) = 2eeP * HSQRT(F), 1/2 =< F< 1 


If K is odd, X = 2*(2P41) * F = 208(2P42) © (F/2), 
HSQRT(X) = 2e8(P4+1) ® HSQRT(F/2), 1/4 =< F/2 < 1/2. 


Let F* = AtF + B, 


103 
104 
! 5 
88 
109 
110 
1 
11 : 
114 ; 
12 3 
HE 
118 : 
120 ; 
121 ; 
158 
124 ; A = 0.453730314(octal), 
125 ; B = 0.327226214(octal), for K even. 
1 § 3 = Av(F/2) + B, 
00 Is 3 A = 0.650117146(octalL), 
BS : 8 3 4 B = 0.230170444(octal), for K odd. 
Py an 
000 189 3 K' =P, for K even 
4 131 : =P +1 for K odd. 
000 133 ; Let YO = 2e*k' * F* as a straight Line spores ination within the 
000 134 ; given interval using coefficients A and B which minimize the 
it : 2 ; absolute error at the midpoint and endpoints. 
B28 130 ; Starting with YO, five Newton-Raphson iterations are performed. 
909 139 ; YEne1] = (1/2) * ( YEn) + X/YEn)) 
909 141 ; The relative error is < 10**-17, 
909 148 ; CALLING SEQUENCE: 
o008 145 : hsart.wh.v = MTHSHSQRT(x.rh.r) 
B83 1h a 
8 143 : CALL MTHSHSQRT(hsqrt.wh.r, x.rh.r) 
00 151: Because an H-floating result cannot be expressed in 64 bits, it is 
00 13 3 returned as the firs orqunent . with the input poreneter displaced 
00 153 ; to the second argument, in accordance with the Procedure Calling 
4 13 : Standard. 
136 : INPUT PARAMETERS: 
aig 138 LONG = 4 : Define Longworg multiplier 
0000008 15 x = 2 * LONG ; Contents of x is the argument 
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MTHSHSQRT ; H Floating Point Square Root routine 16-SEP-1984 01:40:3 AX/VMS Macro V04-00 Page 5 MT 
1-001 MTHSHSORT : Standard H-F Loating SQRT ae 7 91:39:28 MTHRTL.SRC JMTHHSQRT .MAR; 1 . (4) 1- 
0 160 
8 19] :; IMPLICIT INPUTS: None 
Oho 168 3; OUTPUT PARAMETERS: 
00000004 009 192 hsqrt = 1 * LONG ; Contents of hsaqrt is the result 
000 189 3 IMPLICIT OUTPUTS: None 
B68 168 ; 
it 198 3; COMPLETION CODES: None 
00 170; 
4 Ve 3; SIDE EFFECTS: 
0000 158 3 Signais: MTH$_SQUROONEG if x < 0.0 with reserved operand in RO/R3 
0000 174 ; (RO/R1 copied to the signal mechanism vector CHFS$L MCH_RO/R1 by LIBSSIGNAL). 
0000 175 ; Associated message is: "SQUARE ROOT OF NEGATIVE VALUE’ Result is reserved 
44 178 ; operand -0.0 unless a user supplied (or any) error handier changes CHFSL_MCH_RO/R1 
0000 178 ; NOTE: This procedure disables floating point underflow, enables integer 
0000 179 ; overflow, causes no floating overflow or other arithmetic traps, and 
0000 180 3 preserves enables across the call. 
0000 181 ; 
0000 136 ;--- 
0000 18 
0000 184 
41FC 0000 185 eENTRY MTHSHSQRT, ACMASK ; Standard call-by-reference entry 
0002 186 ; Disable DV (and FU), enabie IV 
4 187 MTHSFLAG_JACKET ; Flag that this is a jacket procedure in 
6D 00000000'GF 9E $003 MOVAB G*MTHSSJACKET_HND, (FP) 
0009 ; set handler address to jacket 
0009 ; handler 
0009 ¥ 
0009 188 3; case of an error in special routine 
0009 189 MOVH @x(AP), RO 3; Fetch x from second argument 
OO0E 190 BSBB MTHSHSQRT_R8 : Call kernel HSQRT routine 
0010 191 MOVO RO, ahsart (AP) ; Store result in first argument 
0015 192 RET ; Return to caller 
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MTHSHSQRT ;_H Floating Point Square Root routine 16=-SEP- :3 AX/VMS M v04-00 Page 
1-001 ATHSHSORT_RE - Special HSQRT Routine 6-SEP=1 1984 3 EMTHRIL. SR SRE SMTH HHSQRT.MAR; 1 . 
1 136 -SBTTL MTHSHSQRT_RB - a HSQRT Routine 
16 138 MTHSHSQRT_RB:: 
54 50 70FD 16 19 MOVH RO, R4 ; Test value of x and save in R4/R7 
95 14 OO1A 198 BGTR POs : If x > proceed with computation 
009 31 git 1%) BRW ZERO_NEG : Otherwise handle separately 
O1F 64 3 
O1F 8 3; Here if x > 0. 
bie 50a * 
22 50 =€9 ah 82 POS: BLBC RO, EVEN ; Branch on odd/even exponent 
58 50 eee 8F A3 00 ; 09 ODD: SUBW3 pat 44g RO, R8 ; Save scale factor in R8 
50 FFF 8F BO 00 08 MOVW “X3FFF, RO : Scale x 
54 50 BO 002d 09 MOVW RO. RG : Scale copy of x in R4 
50 50 F6FD 0030 10 CVTHF «RO : Convert x to single precision 
50 13004054 8F 44 0034 11 MULF 2 #LF_ODD_A, RO : cenpute initial approximation 
50 3C4A3F98 8F 40 0038 \¢ ADDF2 #LF ODD-B, : 
20.11 ekg F BRB ITER : Go iterate 
58 50 4000 8F <A3 0044 315 EVEN: SUBW3 #*x4000, RO, R8 ; Save scale factor in R8& 
50 4000 8F B80 O04A 16 MOVW #*x4000, RO 3 Scale x 
54 50 B60 OO4F 217 MOV RQ. oa 3; Scale copy of x in R4& 
50 50 F6FD 0052 $18 CVTHF ; Convert x to single precision 
50 F61A4015 8F 44 0056 19 MULF 2 ag EVEN. ae RO : oneiel first approximation 
50 4B8233FD7 8F 40 8920 $59 ADDF2 #LF_EVEN_B, RO 3 
0064 234 ITERATE: 
0064 2 
0064 se : : : 
0064 $2 3 first iteration - single precision 
0064 § 6 ; 
0064 27 
51 54 F6OFD 0064 228 CVTHF RG, R1 3; RI = 
52 51 50 47 0068 229 DIVF3 RO. R1, R2 3; R2 = 
50 52 40 006C 30 ADDF2 R ; RO = YO + x/Y0 
50 0080 8F A2 O06F 3! SUBW2 #*x80, RO 3; RO = Y1 = (1/2)(Y0 + X/Y0) 
0074 $ ; No overflow possible 
0074 3 
0074 34 ; ' . "= 
0074 35 ; second iteration - single precision 
th $$ : 
- a. oe bond 33 DIVF3 RO, R1, R2 3; R2 = X/Y1 
50 52 40 0078 9 ADDF R2, RO 3 RO = Y1 + X/Y1 
50 0080 8F 4 rt SUBW #*x80, RO 3 RO = Y2 = (1/2)¢¥1 + X/Y1) 
080 6 3 a 
0 8 tz 3; third iteration - double precision 
0 265° 
24 54 76FD 0080 rt} CVTHG R4, R 3 R2/R3 = X 
5 50 99FD 4 4 CVTF RO, R ; RO/R1 = Y2 
52 2 46FD 8 48 DIVG RO, R ; R2 = X/Y¥2 
50 40FD f 49 ADDG R2, R ; RQ = ‘¢ + MNS 
i Ae 50 SUBW #*x10, RO : RO = YS = (1/2)¢¥2 + X/Y¥2) 
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1-001 ATHSHSORT_RE = Special HSORT Routine B-SEP- 1984 $1:55:48 EMTHRTL. § JMTHHSQRT.MAR;1 © (5) 
} | 
: fourth iteration - quadruple precision | 
2 
§ CVTGH RO, RO 3 RO/R3 = 3 
5 DIVH3 RO, R4, -(SP) 3 Stack = x/¥3 
58 ADDH2 (SP), RO : RO = ¥3 + X/Y3 
23 DECW RO : RO = Y4 = (1/2)¢(¥3 + X/Y¥3) | 
61; 
66 ; fifth iteration - quadruple precision 
64 * 
65 DIVH3 RO, R4, (SP) ; Stack = X/Y4 
66 ADDH2 (SP), RO > RO = Y4 + X/Y¥4 
of DECW RO 3 RO = YS = (1/2)¢¥4 + X/Y¥4) 
3 ; scale result by multiplying by scale_factor/2 
71 DIVW : - RE : Divide scale factor by 2 
ie ADDW ; Scale the ag 
7 ADDL Roe Rep : Reset stack point 
a SQRTX: RSB ; Return with feaukt in RO/R3 
6 3; Here if x =< 0 
278 ZERO_NEG: 
79 BEQL SQRTX 3; Return with result = 0 
280 PUSHL (SP) : 
281 MOVZBL #MTHSK_SQUROONEG, -(SP) ; Condition value 
see ASHQ #15, #T, RO ; RO/R3 = result = reserved operand -0.0 
83 ; RO/R1 goes to signal mechanism yector 
see : (CHFSL_MCH_RO/R1) so error handler 
85 3 can modify the result. 
286 CLRQ R2 $ 
287 CALLS #2, G*MTHS$$SIGNAL r Signal error and use real user's PC 
288 3 ependent of CALL vs JSB 
334 RSB 3 eee - RO fockeraa Sean CHF SL_MCH_RO/R1 
291 


Return PC from JSB routine 
| 
| 
j 
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Symbol table 


v04-00 p 
H oan 


acro 
» SRCIMTHHSQRT.MAR; 1 


ACMASK = 4 S1FC 
even 0 G rv 01 

= 
ITERATE $09 64 01 
LF_EVEN_A = F61A4 3 
Hse: SHER 

= 
tone = 96000004 
MTHSSJACKET_HND etereeee xX 4 
RIMEESEGIA™Goosanap act $8 
MTHSHSQRT_RB titty 16 RG 8; 
Sap MOMS sccooaes wa * Bf 
POS QOOOOOI1F R 01 
my i 

= 
ZERO_NEG 00000087 R 01 

¢ eres es eas ee ee wee + 
: Psect synopsis ! 
PSECT name Allocation PSECT No. Attributes 
s OS « 00000000 <( 0.) 00 ¢ O.) NOPIC USR CON ABS LCL NOSHR NOEXE NORD NOWRT NOVEC BYTE 
_MTHSCODE 000000CD ¢«¢ 205.) O1¢ 1.) PIC USR CON’ REL LCL SHR EXE RD NOWRT NOVEC LONG 
Poem renner ere ewe ancane- ~ $ 
; Performance indicators ! 
Phase Page faults CPU Time Elapsed Time 
Initialization 29 00:00:00.07 00:00:01.19 
Command processing 119 Op 8 88 -Be 00:00:04.17 
soetel tebtio sort 8 00:00:00:00 00:00:00:91 
able sor :00:00. :00:00. 
C7 paar. RS BE a 
e ou u $ : e $ : ° 

Psect synopsis otal 3 9 Sha 8 4 00:00:00.02 
Cross-reference output 4% 0.00 ae Se ee 
Assembler run totals 306 00:00:02.26 0:00:12.11 


The_ working set Limit was 900 pages. 

4673 bytes (10 pages) of virtual memory were used to buffer the intermediate code. 

There were 10 pages of symbol table space allocated to hold 19 non-local and 0 local symbols. 
351 source Lines were read in Pass 1, producing 11 object records in Pass 2. 

1 page of virtual memory was used to define 1 macro. 
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! Macro Library statistics ! 


Foor oooooosossessesossesess} 
Macro Library name Macros defined 
_S255SDUA28: (SYS! IBJSTARLET.MLB; 2 SENET 
0 GETS were required to define 0 macros. 
There were no errors, warnings or information messages. 


MACRO/ENABLE=SUPPRESSION/DISABLE=(GLOBAL , TRACEBACK) /LIS=LIS$:MTHHSQRT/OBJ=OBJ$:MTHHSQRT MSRC$:MTHJACKET/UPDATE=(ENHS:MTHJACKET) +MSRC 
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